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' Abstract. Preheating can convert super horizon fluctuations of light scalar fields 

' present at the end of inflation into observable density perturbations. We show in 

detail how lattice field theory simulations and the separate universes approximation 
can be used to calculate these perturbations and make predictions for the nonlinearity 



(N 



D ■ 

. parameter /nl- We also present a simple approximation scheme that can reproduce 

these results analytically. Applying these methods to the massless preheating model, 
• we determine the parameter values that are ruled out by too high levels of non- 

Gaussianity. 

Q ! PACS numbers: 98.80.Cq, 11.15.Kc 
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1. Introduction 
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The non-Gaussianity of primordial density perturbations is a promising way to probe 



lO ! the very early universe. It can be used to determine the precise physics of inflation |I] 

^ I and to distinguish between inflation and its alternatives [21 [3] . Although observations 

I are currently compatible with Gaussian perturbations their sensitivity will improve 

^ I appreciably over the next few years. 

00 i Whilst Gaussian perturbations are characterised completely by their two-point 

^ ' function, it is impossible to give a complete characterisation of non-Gaussian 
perturbations. Non-Gaussianity is usually parameterised phenomenologically by the 

^ I nonlinearity parameter /nl, originally defined [3] for a specific 'local' type of non- 

I Gaussianity by 

C = Co - ^/nl (Co' - (Co')) , (1) 

where C is the curvature perturbation and Co is a Gaussian random field. The definition 
of /nl has since been extended to cover arbitrary non- Gaussian fields ^ [7j, including 
cases in which it is scale dependent. 

Current observational limits from the WMAP 5-year data show that the 
perturbations are close to Gaussian with no departure from Gaussianity [1]: —9 < 
/nl*^^ < 111, —151 < /nl"'^' < 253 at 95% confidence limit. However, a recent work 
[H] using the WMAP three-year data found a significant non-Gaussianity: 26.9 < 
/nl*^^ < 146.7 at 95% confidence limit. Measurements of large scale structure give 
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—29 < /nl'^^ < 69 at 95% confidence limit [9]. These ranges will be significantly reduced 
in the next few years with further WMAP data and data from the Planck satellite. 

To make use of a future observation of non-Gaussianity, reliable predictions must be 
made by each inflationary model. Single field, slow-roll inflation predicts near-Gaussian 
primordial density perturbations, with /nl of the order of the slow-roll parameters e 
and rj pLj. The perturbations generated during slow roll in multi-field models can be 
more non- Gaussian [TUl [TT], [T^ [TB]. 

Even larger values are possible in models with non-equilibrium dynamics during 
reheating f!M [T^ . the process which transfers energy from the infiaton field to matter 
and radiation after inflation has ended. The simplest example of this is the curvaton 
model [ini [IZl [IHl [IS], in which curvaton particles dominate the energy density of the 
universe and therefore influence the expansion of the universe until they decay. 

In many models reheating starts with a brief period of non-equilibrium dynamics 
know as preheating [SUl EI], which results in the rapid production of particles and 
very large occupation numbers. There have been many attempts to calculate the 
perturbations generated using linearised or semilinearised methods in, for example, the 
hybrid inflation model [22], [23] and for chaotic inflation [211 [25], [26] [27] and the separate 
universe approximation [281 [SSI [30] . 

In a recent paper [21] we developed a fully nonlinear calculation using lattice field 
theory simulations and the separate universe approximation. We showed that in the 
massless preheating model [32] parametric resonance can produce very high levels of 
non-Gaussianity, ruling out parameter values. Conversely, other parameter values led 
to strong resonance but negligible non-Gaussianity. 

In this paper we present a more comprehensive study of this model and full details 
of the simulations. 

2. Massless Preheating 

We study a simple variant of chaotic inflation known as massless preheating, the 
dynamics of which have been thoroughly studied previously [2H [251 [32] • The model 
consists of an infiaton field coupled to a massless scalar field x, with the potential 

V{<P,X) = \X<P' + 19W. (2) 

During inflation x is approximately zero and the model behaves the same way as 
the standard single field 0^ chaotic inflation model. We assume that the dominant 
contribution to the density perturbations is generated in the usual way by the infiaton 
field, so that the constraints arising from the linear perturbations are the same as for 
chaotic inflation. This fixes A to be 7 x 10~^^ [33], which we will use throughout this 
paper. This model is actually only marginally compatible with current observations |4|, 
but we still choose to use it because of its convenience: The field dynamics are conformal 
|32j . which means that the relevant physics is at roughly the same comoving scale 
throughout our simulations. As an example, in the case of the potential V{(j), x) = 
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^mcf)^ + ^g^4>'^x^ the inverse mass 1/m sets a fixed (not comoving) physical length scale, 
which has to fit inside the simulation box at all times, requiring much larger lattices than 
those used here. Also, an analysis of the model using second order perturbation 
theory [27] found /nl ^ 0(1000) in the parameter range 1 < g'^/X < 3. 

Unless the coupling ratio g'^/ A is large, the masses of the two fields, = v^3A0 
and = g(f), are comparable during infiation. Furthermore, x is light relative to the 
Hubble rate H except near the end of infiation. Quantum fiuctuations of the x fi^ld 
are amplified and stretched by infiation in the same way as those of the infiaton field. 
At the end of infiation, the x SgIcI will therefore have an approximately scale invariant 
spectrum of perturbations (see Appendix A ). 

The slow roll conditions fail and infiation ends when (p ^ 2^/3Mpl, where 
Mpi = (87rG)-^/2 jg ^j^g reduced Planck mass. After this, the infiaton fields starts 
to oscillate around zero with an amplitude decreasing from this value. The infiaton is 
massless leading the expansion of the universe to be approximately similar to radiation 
domination, a{t) oc t^^"^. During these oscillations, the x field resonates with the infiaton 
field transferring energy away from it, until the amplitude of x becomes so high that 
the dynamics becomes nonlinear [Sni EI]- This process washes out the perturbations 
of X produced during infiation and therefore no isocurvature modes survive. However, 
the perturbations affect the time the resonance lasts and, consequently, the amount of 
expansion during this period of preheating. This means that the perturbations of x 
will leave an imprint in the curvature perturbation, and this contribution is generally 
non-Gaussian. 



3. Separate Universe Approximation 

Our approach [3l] combines classical lattice field theory methods [311 EH] with the widely 
used separate universe approximation [121 ESI EZl EHl EH] ■ It states that points in space 
separated by more than a Hubble distance cannot interact and will therefore evolve 
independently of each other. As long as each Hubble volume is approximately isotropic 
and homogeneous, one can approximate them by separate Friedmann-Robertson- Walker 
(FRW) universes. Gravitational effects are therefore described solely by Friedmann 
equation, which is nonlinear but does not take into account gradients and is therefore 
only valid at long distances. 

In earlier applications of the this approximation [2H1 1221 ED] each separate universe 
was treated as if it was point-like. This means that not only the metric but also the fields 
were assumed to be homogeneous inside each separate universe. As we showed [21] , this 
assumption is generally not valid during preheating. Instead, we allow the fields and 
X to be inhomogeneous on small scales but we will assume that the metric has the usual 
homogeneous and isotropic FRW form. This approximation should be reasonably good 
as long as the size of our separate universes is not much larger than the horizon size, 
but ultimately it should be improved by including short- distance metric perturbations 
in the calculation. It is convenient to choose the separate universes to have a fixed 
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comoving size L. For the FRW metric to be a good approximation for the whole lattice, 
this should be less than the comoving horizon size, L < 1/aH, throughout the whole 
simulation. 

The curvature perturbation ( is given by the perturbation of the logarithm of the 
scale factor a on a constant energy density hypersurface [37J 

C = AN = A\na\H, (3) 

which one can find by solving the Friedmann equation independently for each separate 
universe. If different separate universes have different initial conditions, they evolve 
differently, and this will create curvature perturbations. 

In this paper our emphasis is on non-Gaussianity produced at the end of inflation, 
and therefore we calculate the perturbations during slow roll inflation at linearized level. 
We switch to the separate universe picture when (p = 0ini = 5Mpi, which is a few e- 
foldings before the end of inflation. Using a different value for (p^^i should not change 
our overall results as long as it is in the slow-roll regime. The initial conditions for 
our separate universe calculation are therefore provided by the usual Gaussian field 
perturbations produced during inflation. In a given separate universe, we can write the 
initial field configurations as 

= 0ini + 

= Xini + (4) 
where the mean fields 0ini, Xini are homogeneous over the each individual separate 
universe. They are determined by field perturbations with wavelength longer than L. 
We assume they vary from one separate universe to another with a Gaussian probability 
distribution determined by the power spectra of the fields. The fluctuations 50 and Sx 
correspond to field perturbations with wavelength shorter than L, and therefore they 
are inhomogeneous within a single separate universe. They are represented by Gaussian 
random fields with zero mean and we assume they have the same statistical properties 
in each separate universe. 

As (p aiid X are the only light fields present at the end of inflation, all that needs 
to be considered for the calculation of the curvature perturbation ( is the dependence 
of on their initial values. The inhomogeneous modes 6(f){x) and Sx{x) do not give 
a direct contribution because they have the same statistics in each separate universe. 
Therefore, the curvature perturbation is determined by the initial mean values 0ini and 
Xini- What we therefore need to find is how depends on these two numbers 0ini and 

Xini- 

Our model is symmetric under x ~^ ~Xj and A^ will also possess this symmetry. 
The Taylor expansion of ([3]) can therefore only have even powers of Xim- 

1 d'^N 

C(0ini, Xini) = AAr(0i,i, 0) + -Jf^xli + O(xfni), (S) 

^ "^Xini 

I In the literature it is more common to use lower case delta for this purpose, i.e., 6N. However, we 
reserve S for small-scale perturbations within one single separate universe and use A for large-scale 
perturbations between separate universes. 
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where the second derivative is calculated at Xmi = 0. The first term A^(0ini, 0) is 
independent of Xmi and is therefore exactly the same as in the single-field 0^ model 
where curvature perturbations are known to be highly Gaussian. We will therefore focus 
on the Xini term, which gives a non-Gaussian contribution. To measure the level of this 
non-Gaussianity, we need to find its coefficient d'^N/dXi^i^ which we do by measuring 
the dependence of = In a on Xmi- 

The contribution by the x fi^ld to the non-Gaussianity is not of the simple 'local' 
type (lH), but defining /nl as a suitable ratio of the three-point and two-point correlation 
functions of (, a formula for the effective value of /nl can be derived [3, SHI SI]. Following 
Boubekeur and Lyth [7j the bispectrum can be defined as [12], 



Be ih, k2, ks) = -^/nl [Pc (^i) (^2) + cyclic] . (6) 

We assume that the contribution from the infiaton is practically Gaussian and dominates 
the power spectrum: 

n(^) = n. + nx^n.> (7) 
Be{k,,k2,ks)-Be^{k,,k2,k^). (8) 

It was shown [7] that if the average of x over our currently observable universe is 
negligible and perturbations of both fields are scale invariant, the non-Gaussianity 
parameter is 

For 0^ infiation we would use in this approximation = {H'^/An'^) and = 
(V/247r^eMp[) where e is the slow roll parameter e = 8Mpj/0^. This leads to 

/NL--H^f|^yAM^,ln4r- (10) 

The logarithm refiects an infrared divergence, which is cut off by the length scale at 
which the averages in (Q are computed, i.e., the maximum observable scale. It is a 
result of the assumption that the x field has a flat power spectrum. 

However, the perturbations are not exactly scale invariant, and this changes the 
result. The relevant power spectra are those at the beginning of our simulations, a few 
e-foldings before the end of inflation. It is shown in Appendix A that 

3 , / „ \3(2-g^/X) 



f ^, ^^^P'te) " . .f9VA<(3/2)7V,„. 

g:) AM|,(fi^) (^) e-="^ if,VA>(3/2)Ar_ 



where Nk ~ 60 is the number of e-foldings before the end of inflation when the largest 
scales we observe today left the horizon, and Nsim is the number of e-foldings before the 
end of inflation when we begin our simulations. For our choice of (pi^i this is A^sim = 
Note that we drop the logarithms and other factors of order 1. 

In (P|)- (fTTI) we assumed that the probability distribution of Xini has zero mean, so 
that Xini is zero on average. However, we can only measure density perturbations in our 




Figure 1. The black points are simulation data and the solid red (medium grey) line is 
the order of magnitude estimate given by (j39p for the box size of our simulations. The 
analytic prediction for an infinite box (for which the separate universes approximation 
is, of course, not valid) is given by the dashed red line. Under the assumptions ([7]) and 
dl]) points above the blue (or dark grey) /nl limits (derived from pTjl ) are excluded by 
observations. However, for points above the green (or light grey) lines (derived from 
(|15p ) the same assumptions break down as preheating makes the dominant contribution 
to the power spectrum. This in itself is an interesting and non-trivial possibility. Note 
that at small this limit depends on the total amount of inflation. 



currently observable universe, and perturbations are therefore measured relative to the 
average value xmi over this volume. In order for ( fTTj) to be valid, this average Xini has 
to be small enough. To account for this, we write Xini = X^i + ^Xmi- Then (E]) becomes 



C = Co + C (Xi^ + AXini)' , (12) 

Expanding ( |T2l) and dropping the constant term gives, 

C = 2cx^AXini + cAxL- (14) 

The first term in (|T4l) gives a Gaussian contribution to the curvature perturbation. If 
this has a smaller amplitude than the contribution from the inflaton field, 

2cx^AXini < 10-^ (15) 

where Axini refers to the typical (root mean squared) perturbation at the horizon scale, 
it does not affect observations. Therefore the average Xini can be safely ignored. 



where 
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On the other hand, if the constraint (1151) is not satisfied, preheating is the dominant 
source of curvature perturbations. This means that we would have to use a smaller 
coupling A to obtain the observed amplitude, and other observables such as the spectral 
index would also be different from the usual predictions. 

The typical values of Xmi and Axini are given by the variances 

VAk)ir. (16) 



k 



and 



H 



aoHo 



dk 



(17) 



where Hq is the Hubble parameter today, is the scale factor today and H is the 
Hubble parameter at the start of our simulations. For g^/X < 2 the spectrum, (k), 
is red tilted (there is more power at larger scales) and the integral is infrared divergent. 
As we show in Appendix A[ this means that the constraint depends on the total amount 
of inflation A^tot- For g'^/X — > 0, we flnd 



97r2 
4A 



{N,,,No)-^/' X 10-5 ^ 3000 X {N,,,/100)~'/\ 



where Nq ^ 60 is the number of e-foldings of inflation after largest currently observables 
scales left the horizon. For g'^/X > 2 the spectrum is blue tilted (there is less power at 
larger scales) and the integral converges. Because of this, the constraint on c becomes 
rapidly much weaker, as shown in flgure [H See Appendix A for further discussion of 
this issue. 



4. Analytic Approximation 

In this section we derive an analytic approximation of the second derivative d'^N/dXi^i, 
and therefore /nl, by linearising the fleld dynamics of preheating. This will give an 
order of magnitude estimate to compare with the fleld theory simulations described in 
section 

The dynamics of the flelds are described by a system of two coupled ODEs: 

4> + 3H<p- + + ^Vx' = 0, (19) 

X + 3Hx-\v\ + g'(j)\ =0. (20) 

The details of the fleld dynamics are given in detail in [32], but we will repeat some of 
the more important points here. At the end of inflation the inflaton fleld, (p, reaches 
the bottom of its potential and starts to oscillate with decreasing amplitude. During 
and immediately after inflation the x fi^ld is approximately zero and we can ignore the 
terms involving x in (IT^ . along with the gradient term as the fleld is approximately 
homogeneous. In terms of rescaled fleld (p = acf) and rescaled conformal time f deflned 
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Figure 2. Resonance structure of the massless preheating model. The contours show 
the Floquet index, /i, from the solution of (|22p . The peak of the first band is at 
/\ — 1.875. The four marked sections are shown in figure [3] 



hy dr = a ^X^/'^A^dt^ the oscillations are approximately described by the Jacobi cosine 
functioifl [13], 



0(f) = A<^cn(f,l/V2), (21) 

where is the constant amplitude of the oscillations. This is set to = 2a/3Mpi, 
the value of when slow roll is violated. The oscillations in the field give rise to an 
oscillatory mass term for the x field. At linear level, a Fourier mode of the rescaled field 
X = ax with wave number k satisfies the Lame equation. 



Xk + 



2 T fc2 



.2 , 9__2t^ 

X 



+ ^cn^(r, 1 



Xk = 0, «' = ^- (22) 

In the space of the two constant parameters and g"^/ X, ( l22l) has resonance bands 
in which the solution grows exponentially, 

Ur) = e'^^^'^'/'^Vlr), (23) 

where /i(K, gf^/A) is known as the Floquet index and /(f) is a periodic function. This 
means that energy is transferred rapidly from the infiaton, 0, to the x field. The Floquet 

§ We follow the usual convention in cosmology and some mathematical works |43| and use 1 / ^2 in the 
second argument of the Jacobi cosine function. Note that the same function is often given as cn('r|i) 

m. 
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Figure 3. Cross-sections of tlie resonance structure shown in figure [5] at constant 
values of /\. 



index for a particular set of and / A can be calculated by numerically solving (l22l) 
over one period of oscillation and finding the eigenvalues of the matrix that relates the 
final and initial values. 

As figure [2] shows, the resonance bands stretch through the {g'^/\,H?) plane 
diagonally so that for any value of g"^ /\ there are some resonant modes. Figure [3] shows 
how the Floquet index depends on k at some chosen values of g'^jX. When g"^ <^ A, 
the bands are narrow and the resonance is weak. Effective preheating therefore requires 

gV^ > 1- 

It is important for our purposes here to note when the modes of the x fi^ld grow 
on large scales. This is when modes with = fall in the resonance bands in figure 
121 The edges of these resonance bands fall at g'^/X = |n (n + 1) [32]; the first band is 
1 < g'^/X < 3, the second is 6 < g^/X < 10 and so on. 

To estimate the curvature perturbation produced by the resonance, we describe the 
resonance using the linear approximation. When the amplitude of x has grown so much 
that this approximation is no longer valid, we assume that the resonance stops and the 
fields equilibrate instantaneously. 

In the first instance we assume that the universe expands with the radiation- 
dominated equation of state, so that a oc t^/^. In terms of rescaled conformal time 
f, the scale factor is 



a r 



1 



A, 



-T. 



(24) 



12Mpi 

In the following we will use the scale factor a as the time coordinate rather than f. 
In the linear approximation field perturbations remain Gaussian. This means 
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that the modes are independent and the zero mode and inhomogeneous modes can 
be separated, 

{x') {a) =f {a) + {5x') {a) , (25) 

where x = (^X- 

Following (123|) . we approximate the evolution of mode with comoving wavenumber 
k by Xfc(a) oc e'^'='^, where fx = \/l2fj,Mpi/A^ is a rescaled Floquet index. Consequently, 
the first term in (1251) grows as 

f (a) = f (1) e2'^°('^-^) = x^^.e'^^'^^-'l (26) 

Note that for our choice of ajni = 1, Xini = Xini so we drop the tilde. We approximate 
the evolution of the second term by, 

(Sx^) (a) = j (e2A.{-i) _ 1) ^ ^2^2ii^Ma-i)^ (27) 

where = XA^ and /imax is the Floquet index of the mode with the largest Floquet 
index for a choice of g'^/X. More precisely, the non-exponential prefactor is determined 
by the shape of the resonance band and it therefore depends on the coupling ratio g"^/ X. 
It would even be reasonably straightforward to compute it numerically by evaluating 
the integral. However, since the only relevant dimensionful scale is m, this serves as a 
sufficient order-of-magnitude estimate. 

Also by inspection of (JTOll it can be seen that the system becomes nonlinear at scale 
factor ttni when g'^ (x^(ctni)) — ~ (see also figure [5]). Substituting fl26l) and fl271) 
into (P^j) gives, 

- 9^ (xfrne^^"""' + m^e^'^-^^""') . (28) 

The nonlinearity scale factor depends on Xini, and to determine it we Taylor 
expand it as 

ani(Xini) = (1 + CnixLi)ani(O) + 0{xL)- (29) 
We first solve the equation for Xini = 0, and find the corresponding value ani(O) 

~ ^2^2g2/i.„axa„l(0) ^ ^^^(Q^) _ _]_ 1 (3Q^ 

/^max 9 



We then determine Cni by substituting (1^ to 

~ g^ (x?„ie'^"""'^°^ + m'e'''— "-'("^e^'^— + O(xfni)- (31) 
Rearranging this gives 

rs^ 3^ g2(/io-/imax)anl(0) (22) 

2/imaxanl(0)m2 

and substituting for ani(O) from (130|) and = AA^, 

1 g^ /OQN 

<p a 
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Equations (!29l) and (!33l) tell us the value of the scale factor at the time when 
the dynamics becomes nonlinear. This, however, does not yet give the curvature 
perturbation. According to ([3]), the curvature perturbation is given by the scale factor 
at some fixed energy density or, equivalently, some fixed Hubble rate, which we denote 
by if*. The result does not depend on the choice of H^: as long as it is chosen to be 
after preheating has ended. 

To calculate the curvature perturbation from fpUl) . we assume that once the 
dynamics has become nonlinear, the equation of state is exactly that of radiation, and 
therefore H{a) = {a^Ja^)H^i, where H^i is the Hubble rate at ttni- Inverting this we 
find a{H^) = Hn\/ H^ani- This implies 

A In a(ifj = -Alni^ni + AlnOni = ( 1 + -^t"-" ) AlnOni. (34) 



2 V 2 dlna^ 

The second derivative at constant if in ([5]) is therefore 

d^N f ld\nH\ , , 

If the equation of state during the resonance were exactly that of radiation, the two 
terms inside the brackets would cancel exactly and there would be no contribution to 
the curvature perturbation ([3]). However, the expansion rate depends on the phase of 
the oscillations, and in fact we have 

2 / ~\ 2 



d\wE 1 



d\wa 2M|j \ da 




(36) 



Using (1211) . we find 



— «-6(-cu(f,l/y^)) . (37) 

This oscillates between —3 and 0, and therefore 

1 / ldlnii\ , , 

- - < 1 + — < 1. 38) 

2 - V 2 d\na J - ^ ' 

This demonstrates that the two terms in fl35p do not generally cancel, and that in 
fact the conversion factor is or order one but can have either sign. Therefore we adopt 
as our analytic prediction 

_.±2,„.±_^,-^. (39) 

This is shown in figure [1] and its shape can be understood by comparing with the 
Floquet indices shown in table [TJ The effect is only present when the zero mode is 
within a resonance band, i.e., /i > 0. Near the lower end of the resonance band, when 
1 < < 1.4 (see the blue long-dashed line in figure [3]), the zero mode resonates but 
the maximum Floquet index is achieved at k > 0, and therefore /io//imax < 1- 

Above this, for 1.4 < g'^/X < 2.986 the resonance is at its strongest at k = 0. 
For the solid red line in figure [T] we have chosen to use as /Umax the value for the 
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longest wavelength available in our finite simulation box, to allow direct comparison 
with our numerical results. However, this makes the result dependent on the box size 
and therefore obviously unphysical. Taking the naive infinite- volume limit L oo would 
give us /io//^max = 1 (fed dashed line), but this would violate the constraint L < H^^. 
Therefore, we have to accept that our approach is not reliable in this range of g'^/X 
and that the true value is probably somewhere between these two limits. Nevertheless, 
setting /io//imax = 1 gives a useful rule of thumb: d'^N/ (9xfni ~ 1/A when g"^ / X is close to 
the top of a resonance band. Therefore, we can see that d'^N/ c^xf^; will be large because 
infiation constrains A to be small (~ 10"^"^). 

Finally, in the range g'^/X > 2.986, the strongest resonance is found in the second 
resonance band (pink dot-dashed line in figure [3l), and therefore the prediction falls very 
steeply in this range. 

5. Simulations 

In the simulations we employ three-dimensional classical field theory lattices. This is a 
standard method of solving such systems [HH [35], US, HHl HTJ HH], but before [31], it had 
not been used in the context of the separate universe approximation. To determine the 
expansion of the universe, we couple the equations of motion (|T^ with the Friedmann 
equation 
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Figure 5. Evolution of the x field during one simulation for /X — 2.7. 



where the energy density is calculated as the average energy density in the simulation 
box, 



P 



1 

Z3 



1 

12* 



1 -2 1 



(41) 



The coupled system of equations (fT^ and (HUj) are solved on a comoving lattice with 
periodic boundary conditions. We solved them in conformal time {dr = a~^dt) using 
the second-order Runge-Kutta algorithm for the field equations coupled to an Euler 
method for the Friedmann equation. (Details of the algorithm used are presented 
Appendix Bj). After each Runge-Kutta timestep for the fields the integral (HI 



m 



IS 



performed and the Euler timestep for a is made. While it is possible to use a Runge- 
Kutta system to solve all the equations, we found the different order of the equations 
( ffT^ is second order and fHUl) is first order) leads to a cumulative error in a and, 
therefore, numerical errors in the final results. This is not the case with the Runge- 
Kutta-Euler hybrid algorithm used here. 

The separate universes approximation has been applied to massless preheating 
previously [28], EHl [30,119], but without including the field gradient terms in (fT9l) . In these 
works the initial x is varied and the log of the scale factor, A^, at some later H is found. 
These works found the function A^ (Xini) to be random, suggesting chaotic dynamics. 
However, the calculation does not include contributions from the inhomogeneous modes. 
It is the averaged dynamics which lead the function A^ (Xini) \h to depend on an average 
across all of the modes. As we shall see this addition of additional degrees of freedom 
leads to the chaotic dynamics being smoothed for sufficiently small values of Xim- 

In setting up the initial conditions, we treat the homogeneous and inhomogeneous 
modes differently. The lattice average of x is set equal to Xini- For the inhomogeneous 
modes, we follow the standard approach [Mll35l H5l H61ITF] . The x fisld is given random 
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Figure 6. Evolution of the scale factor during one simulation. We begin the system 
a few (~ 3) e-foldings before the end of inflation. During preheating the scale factor 
evolves approximately as in radiation domination. 



initial conditions from a Gaussian distribution whose two-point functions are the same 
as those in the tree-level quantum vacuum state, 

1 



k 



(XfcX,) = (27rmfc + g)^|, (42) 



where ujk = \f^^^rrn^ = -^/fc^ + g'^(f)f^i- All other two-point correlators vanish. Early in 
the simulation (at the end of inflation and the beginning of preheating) the dynamics 
are linear, and the time evolution of this classical ensemble is identical to that in the 
quantum theory. Later, when the field evolution becomes nonlinear the occupation 
numbers are so large that the classical theory is also a good approximation to the 
quantum field dynamics. 

The inhomogeneous modes of are populated similarly to x- The initial conditions 
for homogeneous mode of the field is 0ini = 5Mpi. This is sufficient to drive ~ 3 
e-foldings of inflation. The initial scale factor is a = 1. The usual slow roll techniques 
(which are valid only during this early part of the simulation where (p > 2v^Mpi) are 
then employed to give the initial flrst derivative of with respect to conformal time r 
to be ((i0/(ir)ini = — ^AVMpi0ini. From this point the Runge-Kutta-Euler algorithm 
described in Appendix B is used to solve (fT9l) and (HOl) . 

FigureHlshows the evolution of the mean of the inflaton fleld during one simulation 
time. The fleld begins above Mpi which drives the initial inflation which can be seen 
in flgure El It rolls to the bottom of the potential and proceeds to oscillate about the 
minimum. Example output for the x fi^ld is shown in flgure 

We flxed the inflaton self coupling to A = 7 x 10~^^, which, assuming that 
inflation makes the dominant contribution to the power spectrum, gives the observed 
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Figure 7. The evolution of the log of the scale factor, ln(a), at the end of a sample 
of simulations for / X = 1.875, with radiation-domination-like expansion divided out. 
Each line averaged over 100 runs. The three curves arc for different Xini- Dashed: 
Xini = 0. Dot-dashed: Xini = 1.0 x 10^'^. Solid: Xini = 1-9 x 10^^. At late times these 
tend back towards radiation domination (a horizontal line in this plot). In the inset it 
can be seen that the lines a parallel at late times, resulting in the perturbations being 
locked in. 

level of CMB fluctuations [I33j. We used lattices of 32^ points with comoving spacing 
Sx = 1.25 X 10^ and time step 6t = 4x 10^ in Planck units. The lattice size is smaller than 
in most other studies of preheating because the calculation of curvature perturbation 
requires a large number of runs for each set of parameters, whereas in other, much less 
ambitious studies a few runs have been sufficient. It should also be noted that the 
constraint L < 1/aH limits the size of our lattice, and even with the current choice, 
the comoving horizon size is briefly somewhat smaller than L. Therefore, the number of 
points could only be increased by reducing 6x, i.e., making the lattice flner. However, 
the relevant physical scales are longer than 6x, and we see no indication that this would 
improve our results. 

It should be noted that the gradient energy of the 'quantum' fluctuations f H2l) gives 
an ultraviolet vacuum energy contribution to the energy density. This can be estimated 
to be puv ~ Sx~'^ ~ 10~^^, and we have to make sure that this is much less than the 
physical energy density stored in the inflaton potential Pphys ~ i-^^fni ~ 10~^^. 

According to ([5]), the contribution to the curvature perturbation from preheating 
depends on the second derivative d'^N/dxini constant H = H^, where Xini is the value 
of the zero mode at the beginning of the simulation. If we calculate the second derivative 
at late enough times, when the system has reached a quasi-equilibrium state, the result 
should be independent of because the curvature perturbation is conserved. 

Our simulations give us the functions a(r) and H{t), from which we obtain a{H). 
This is shown in flgure[7l where we have subtracted the underlying radiation-dominated 
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Figure 8. The dependence of N on Xini for g^/X = 1.875 measured at H = 
5.53 X 10^^^. The curve shows a quadratic fit for low N on Xini to the function 
-^(Xini) = -^(0) + cXhii- Due to the symmetry of the system simulations are only run 
for positive Xini for all other /X. 

evolution a oc H^^/'^. Initially, we can see the effect of the coherent oscillations of the (f) 
field, which correspond to (1371) . When the dynamics become nonlinear at H 10~^°Mpi, 
these die away. However, as the inset shows, smaller oscillations due to non-equilibrium 
effects and statistical errors remain. 

To find d'^N/dXini we need the scale factor for different Xini at some chosen value 
H = H^:, whereas the output from the code is at discrete values of H which are different 
for each run. Therefore we have to interpolate the data from each run to H = H^,. 
To simultaneously remove the effect of the transient oscillations, we fit a power-law 
function to the data over a range oil/ H that is longer than the characteristic length of 
the transient oscillations, and use that to determine = \\ia{H^) for each Xim- 

From this data, we obtain the second derivative d'^N/dXi^i by doing a fit at low Xini 
with a quadratic function, 

iV(Xini)=iV(0) + cxL, (43) 

so that d'^N/dxfni = 2c. For each Xini we repeated the simulation between 60 and 
240 times, each with a different random realisation of the initial fluctuations. The 
averages for each xim can then be plotted (see for example figure [H]) and a best fit for 
the parameters iV(0) and c in ( H3|) can be found. As the figure shows, the function ( l43l) 
is only a good fit for small Xim- We must therefore make a choice of which points to 
do the fit over. To do this, we start with a small number of points included in the fit 
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Figure 9. The evolution of the fit parameter c in figure [5] and for two other choices of 
g^/X. FiUed circles: g^/X = 2.7. Unfilled circles: g^/X = 1.875. Crosses: g^/X ^ 1.175. 
Note that c does not change at late times. 



and steadily add more points to the fit. When the statistical psr degree of freedom 
grows significantly above one, points of that Xmi and above are not usedlJi] 

In figure we show how the result depends on the value of H^, at which it is 
measured. The plot confirms that at late enough times the result is independent of 
as it should be. 

The final results measured at = 5 x 10~^^Mpi are shown in table [1] and in 
figure [H Table [1] shows that the sign of the result varies as suggested by the analytic 
result (l39l) . The figure also shows the analytic result in ( l39l) . which is in very good 
agreement with our data, particularly in the first resonance band. This demonstrates 
that the calculation in Section H] captures the relevant physics, and gives a potentially 
very useful way of calculating curvature perturbations in other models without having 
to carry out numerical simulations. Because of this, and because the analytic result 
covers all values of g"^/ A, we use it to draw some further conclusions. 

In figure [1] we also show where the amplitude and the non-Gaussianity of the 
perturbations exceed the observed values according to (|T^ and (|TT|) . Comparing the 
data to the non-Gaussianity limits shows that when the large scale modes are within 
most of the first resonance band (1 < g^/X < 3) the prediction for /nl is far outside 
observational bounds. Elsewhere the predicted /nl is negligibly small. The ranges 
of g'^/X where /nl is compatible with current data but potentially observable, i.e., 
1 ^ I/nl| < 100, are extremely narrow: 1.0321 < g^X < 1.0408 and 2.9934 < g^/X < 
2.9941. 



II There is overlapping notation here, 
and Xini refers to the scalar field x- 



'x^ per degree of freedom' refers to the statistical technique 



Non-Gaussianity from massless preheating 



18 



Table 1. The results from simulations are shown in the final column. The 
corresponding heights and positions in the resonance structure (figure [5]) are shown 
for reference, /imax is the largest Floquet index for non-zero k"^. Note that for many 
simulations this is at the largest scale in the lattice. 
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Within most of the first resonance band (1.060 < S'^/A < 2.992) the constraint 
f|T5|) is not satisfied, meaning that the amphtude of the power spectrum due to the 
hnear term in ffT^ also exceeds the observed value. This means that also the Gaussian 
perturbations are dominated by the contribution from preheating. By varying A their 
amplitude can be tuned to the observed level, and we plan to investigate this interesting 
possibility further in a future work. However, in this paper we simply conclude that for 
the usual choice of A these values of g^/X are ruled out. 

Figure [T] also shows that when the zero mode falls in the second resonance band, 
6 < g"^/ A < 10, and higher resonace bands, the constraint (fTSll is satisfied but preheating 
does not lead to observable non-Gaussianity. 

As discussed in section HI our method has its limitations: the sharp spike close to 
g'^/X = 3, which is present both in the analytic and numerical results, is unphysical and 
can be interpreted as a finite-size effect. As shown in figures [2] and [3l for these values of 
g'^/X the first peak in the Floquet index /i is narrow and includes only a few modes with 
low K. It is these modes which make the dominant contribution to d'^N/dxl^^, and as 
g'^ /X approaches 3 our finite simulation box contains fewer and fewer modes which are 
in resonance until the zero mode is artificially isolated by being the only mode falling 
in the first resonance band. This problem cannot be solved by using larger lattices, 
because if the lattice size L exceeds the horizon size we cannot assume that the 

whole lattice is described by a homogeneous FRW metric and the results ( ITOl) and (ITTll 
would not hold [3T]. The dashed red line in figure [1] shows the prediction from fl39|) for 
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a box of infinite (and therefore unrealistic) size. The closeness of this to the prediction 
for a box of the size of our lattice simulations and to the results of the simulations 
shows that our conclusions do not depend heavily on the box size apart from in these 
sharp spikes. However, problems arising from the box size being slightly larger than the 
horizon could be addressed by introducing linear metric perturbations on the lattice as 
was done in [H], as long as deviations from homogeneity within the lattice are small. 

6. Conclusions 

In this paper we have presented full details of the method for calculating the curvature 
perturbations produced by non-equilibrium dynamics after the end of inflation which 
we introduced in [HT]. The method can be applied to many models, and we have 
demonstrated it by analysing the massless preheating model. 

Our results show that preheating can have a very large effect on the curvature 
perturbations, but that it depends sensitively on the model and its parameters. In the 
case of massless preheating, the contribution is small if g^/X is large, because the x field 
is massive and its fluctuations are suppressed. For g'^/X ~ 0(1) we find a non-trivial and 
interesting structure. If preheating is dominated by long-wavelength modes, it produces 
large non-Gaussian curvature perturbations which are incompatible with observational 
limits on /nl- For most of these values even the amplitude of the perturbations is too 
high. The amplitude can be reduced to an acceptable level by decreasing the value of 
A, leading to a scenario in which the observed perturbations arise predominantly from 
preheating. We will investigate this interesting possibility in a future publication. 

There are, however, narrow regions near g"^/ A = 1 and g"^/ A = 3 that are compatible 
with current observations [1] but for which /nl is large enough to be observed with future 
experiments such as the Planck satellite and can even saturate the current observational 
bounds. 

It is important that preheating, a small addition to the simplest inflationary models 
which was originally introduced for very different reasons, can produce such levels of 
non-Gaussianity when slow roll inflation alone cannot. Unfortunately it also means that 
observation of non-Gaussianity can neither prove nor disprove these models. 

What is also important is that even if x is light, preheating produces no observable 
effects if it is dominated by inhomogeneous modes. These values of g'^/X are allowed in 
spite of the presence of isocurvature x field fiuctuations, because they are wiped out by 
preheating. 

It will be very interesting to see how well our findings generalize to other, more 
realistic inflationary models with preheating or other non-equilibrium phenomena. Our 
numerical method can be readily applied to any bosonic model, although this will be 
somewhat more costly since massless preheating has some special properties that make 
simulations particularly easy. Perhaps an easier way to study such models would be to 
use the analytic approximation we present in section HI which reproduces our numerical 
results and does not require numerical simulations. 
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Appendix A. Power spectra at the beginning of simulations 
Appendix A.l. Amplitude of perturbations 

The curvature perturbation generated during preheating is a function of the average 
value of the field x over the simulation volume at the start of the simulation. The 
spectrum and other statistical properties of the curvature perturbation are therefore 
determined by the spectrum of x at the time when them simulation starts, which is 
shortly before the end of inflation. 

During inflation the evolution of 6xk is given by 

6xk + SH6xk + 9'<p''Sxk = 0. (A.l) 
If we approximate the field to be massless, (^^0^ <^ {3H/2)^ then the modes are frozen 
once they leave the horizon and we have the standard result for a massless field [33], 

^ - " ^ (A.2) 



27r^ 



\SXk 



4vr2 



k=aH 



16 

y 



(A.3) 



and 



(A.4) 

In these equations, measures the number of e-foldings before the end of inflation, 
but it should be noted that inflation actually ends at = 3/2 where the slow roll 
conditions fail as rj = 1. The slow decrease of H in flA.3p makes the spectrum flA.2p 
scale dependent. 

The X field is also not exactly massless during inflation, which introduces another 
source of scale dependence to V^. The Hubble damping term 3H in (lA.ip decreases 
more rapidly during inflation than the mass g(f). At N = Ncru = (2/3) (7^/ A, when 
g(j) = 3H/2, X becomes underdamped. A mode leaving the horizon before this will 
experience Nk — Ncrit e-foldings of overdamped freeze-out followed by A'crit e-foldings 
of underdamped oscillations, where Nj^ is the number of e-foldings before the end of 
inflation when the mode leaves the horizon [50]. Modes leaving the horizon during the 
short underdamped period can be ignored, since they are never amplified by inflation. 

Our simulation starts at A^ = Asim- If Agim > Acrit, which corresponds to small 
g'^/X, then the entire underdamped period is calculated by our lattice simulations. In 
this case, the overdamped period adds a factor of \^iU\ 



(A.5) 
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to dAJ). Substituting, flA^l) and fIXil) . 







exp 









9 3 
4 ~ 2A"iV ~ 2 



and cit 

2 



leads to, 



Solving the integral, ( 1A.2I) becomes [50j . 



where 



= Nk - iV3i 

+A^crit log 



3F(Affc,Ar,in,) 



(A.6) 



(A.7) 



k=aH 



crit 



, , , . (A.8) 

In the case where Ngim < Nait, i-e., for large g'^/X, there is a similar overdamped 
contribution from A^^^ to A'crit and also a contribution from the underdamped period 
between A'crit and A'sim 1501, 



giving 



exp 



^cnd 



3H 



dt 



exp 



N ■ n 2 

dN 



~3F{Nk,N,,it)-2S-+3N,i„ 



(A.9) 



(A.IO) 



k=aH 



As discussed in Section |3l our assumption that the curvature perturbations are 
dominated by the infiaton field leads to the constraint (fTSll on the typical values of Xmi 
and Axini. These are given by (fT6l) and (|T71) . Changing the integration variable from k 
to Nk, we can write the equations as 



and 



No 



1 



No / ^ 



A^. 



ciAT. 



dNk, 



(A.ll) 



(A.12) 



where A^ is the number of e-foldings before the end of inflation, A"o ~ 60 is when the 
largest currently observable scales left the horizon, and the cutoff A^tot > is the total 
number of e-foldings of inflation. Evaluating these integrals numerically, we find the 
constraint shown in figure [H 

At small g^/X the integral diverges as A^tot —>■ oo, and therefore the constraint is 
depends on the total number of e-foldings A'tot- In the limit g'^/X —>■ 0, we have A'crit = 
and F{Nk, A'sim) = 0, so that 

Pxik) = S ^ AaM^^AT^. (A.13) 



47r2 
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We can therefore solve the integrals (1A.11I) and (1A.12P easily and find 

(Axini') ^ ^AMliiVo^ (A.14) 
Substituting these into (|T5i) leads to ( |T8|) . 

Appendix A. 2. Non-Gaussianity 

When making estimations of /nl from © we also need to know the power spectrum 
of the infiaton (j), which we assume to make the dominant contribution to the curvature 
perturbation. Perturbations of (p evolve according to 

6'^k + SH6^k + SX(j)^S(f)k = 0, (A.15) 

which is identical to (lA.ip if we replace 3A. This corresponds to A^crit = 2, and 

therefore the power spectrum is given by the analog of (1A.7I1 . 



tt2 



(A. 16) 

k=aH 



-3F(Affe,Ar,in,)|^ 

47r2' 

In order to estimate /nl we must modify ([9]) which was derived using the method 
described in [7j for the case of non-scale-invariant power spectra. We split this into 
two parts. Firstly the we find the form of the ratio of power spectra V^/V^ and then 
secondly we modify logarithmic factor. 

For the case in which Asim > ^crit (the simulation begins before x fi^ld enters the 
underdamped period) we substitute ( 1A.7I) and ( 1A.16P in to V^/V"^ using the slow roll 
solutions to estimate the scale-invariant component, 

Making the approximations A^^ ^ Acrit and Nk ^ Agim the function F{Nk, A'sim) 
becomes 



F{Nj„N,,^) = AT^itlog I 1 , (A.18) 



and therefore, 



A^sim 



V' 16 . , / N, ^ '^'-''/'^ 



X 



In the case in which A'sim < ^crit (the simulation begins after the x fi^ld has entered 
the underdamped period) we substitute ( lA.lOp and ( 1A.16P in to V^/V^: 

T)^ 1 « 2 



which in the large A^^ limit is 



-^AM^Y^V" 7^ ) e-. (A.21) 



Sim 
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The origin of the logarithm in ([9]) is in the estimation of the integral [7j 

r I ,3 I A |3 ^x(g)^x(g - fci)^x(g + ^2). (A.22) 

^aoHf, 0' |g-fci| |g + fc2| 

We approximate fci ~ A;2 ~ ^- For large g the integral goes as ~ -rj, so we assume that 
the dominant contribution come from g -C /c and the integral becomes 

A>x(^), (A.23) 

which in terms of is 

' Xl'^j / AAJ T> ( AT \ \ ^ ^ 



^ J„/WiV,)^l-^j. (A.24) 

At large we can drop the ^1 — -^-^ term and the power spectrum goes as 



AT,, 



~ / rfAT.AT^^ I _^ J ^ I . _ (A.25) 

From this it can be seen that as found previously [50] the spectrum is only scale invariant 
for g"^/ A ~ 2 in which case (1A.22P is proportional to 

dNk ~ A^o - iVfc ~ In ( ^ j . (A.26) 



In general we have, 

Combining this with and in the A"sim > A'crit case with (1A.19I1 we have 



/nl ^ 




In the A'sim < A'crit case we substitute (1A.21I) into (I9l): 



-(log^V^^^ |. (A.28) 




(A.29) 



Dropping the numerical constant and the logarithms leads to ([TT 



Non-Gaussianity from massless preheating 



24 



Appendix B. Numerical algorithm 

Here we will present the numerical algorithm used in our lattice simulation to evolve the 
field equations and the Friedmann equations. This is done by the numerical integration 
of ( fT9l) and f HOl) . The variables to be solved for are 0*''^(r), x^^^ij) ^iid a(r), for all i, 
i and k where i G {0, S — 1}, j G {0, S' — 1} and k G {0, 5* — 1} indicate the 
position in the lattice, and there are S'^ points in the lattice. We use conformal time 
r defined by dr = dt/a{t), and ' will indicate the derivative with respect to r. The 
conformal time step is denoted by 6t and the comoving lattice spacing by 6x. 

The field and Friedmann equations are coupled in a leapfrog-like fashion, by defining 
the scale factor a at half-way between time steps of the field evolution. To achieve this, 
the scale factor is first evolved by half a timestep using slow roll. Making the usual 

assumptions, a'(rini) = y lot'f' Q-(''"ini)^; and then. 



a ^Tini + = a(rini) + Ya'(rini). (B.l) 

Slow roll is not imposed after this initial half-step. The fields are then evolved one 
timestep to be half a timestep ahead of ^(t), and then a(r) is evolved one timestep to 
be half a timestep ahead of the fields, and so on. 

The field evolution is by a standard fourth order Runge-Kutta method [Bj. The 
first derivatives are given the status of independent variables: 0p^^ = ((/)*-'*^)' and 
Xp = {x^-''') ■ Therefore at the nth timestep: 

= - 2^0f„ - al^, (a ' + / {x'^tY) (B.2) 

(^ijk\" _ (^ijk\' _ r ( jSjk lijk ijk ^ijk \ 

= VY^i - 2-^xfn - <^i9' (0^)' X^, (B.3) 



2 



where for X G {0, x\i 

^2-w-ijk _ ^ / , -w-li-iyk , v-i[j+l]k 

^ n " ' n ' n 

^ ^ x^^'lf-^l - 6X^^'^) . (B.4) 

The square brackets indicate addition or subtraction in modulus S. We the define the 
Runge-Kutta parameters: 

1^1 = 0p't 

i^ = uif'i<i^ix^',x^i:) 

f>ijk _ ijk 

-'^13 ~ Xp n 

ryijk _ r /lijk lijk ijk 'ijk \ 

-n-U —Jx\y^ n'^pn^Ap J Xp n) 

pijk _ ±ijk I 1 pijfc 
-"■21 — 'Tp 2 
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-"-22 — /</>lV^ ll'^P" 2 P ^2^3'Apn+2^4 ) 

1 

yijk — f , ( Aijk I 1 riijk lijk , ^ r>ijk ijk , ^ pijk ijk , ^ pii^ 



2 



i?3^ = X^t + 2^24 

i?3^ - /x ('/•^^n + \Rii^. 'Pfn + ^^2^, xf + ^^^2!^ , Xp't. + ^^^^ 
= + Rf, + ^3^, Xf + Rll Xfn + Rf) 

D«ifc _ ^ijk I pjjfe 

-"-43 ~ Xp n "T -"-34 

^4^' = fx + Rf: ^fn + R!i. xf + i^l', Xf„ + i^^') ■ (B.5) 

Prom these we can find, 

xf„+, = xf„ + y (flS? + 2< + 2B-y + flf) . (B.6) 
The evolution of a is by an Euler method. At each timestep we have the sum, 

i=0 j=0 fc=0 V n-i n-i 

where for X G {0, x}, 

VX^^'^ = (Xl^+i]^*^ + Xl^+i]^*^ + X^I^+i]'^ + X^'V^l'^ + + X^^jf-^]) . (B.8) 

Putting this into the Priedmann equation gives. 
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preheating [JCAP 08 (2008) 002] 
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Our simulation code, which was used to produce the numerical results in Fig. 1 
and Table 1, contained an error. The lattice discretization (B.8) of the gradient term 
in the energy density (41) was incorrect, and the calculated quantity was therefore not 
the conserved energy that corresponds to the discretized lattice field equations (B.2) 
and (B.3). Equation (B.8) also contained unrelated typographical errors. The lattice 
gradient we used in the simulations was 
1 

\jvijk _ ( v[i+l\jk _ v'[i-l]jk Yi[j+l]k _ _ X''i-i\k-^\\ (-\\ 

where 

(2) 

The correct discretisation, consistent with equations (B.2) and (B.3), is 

VX'^l = ^ (Xl^+^l^'^ - X'^'^, Xt^'+i]-''^ - X*^'^, X^^'[f+il - X^^'^) . (3) 

As a consequence of this error the numerical results cannot be trusted. This does 
not affect the validity of the method we described in our paper. Corrected simulation 
results will appear in Ref. [1]. We thank the authors of that paper for pointing out this 
error to us. 
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